################################################################################
## TITLE:		  poll_regulation
## DATA:      ca2019_ces_poll_regulation.dta
## AUTHORS:   Review version
## PAPER:     The Regulation of Pre-Election Polls: A Citizen’s Perspective
## DATE:		  February 2023
################################################################################

getwd()

setwd("C:/Users/Utilisateur/Desktop/ijpor")

## INSTALL AND LOAD PACKAGES

#install.packages("nnet")
#install.packages("effects")
#install.packages("magrittr")
#install.packages("ggplot2")
#install.packages("tidyverse")
#install.packages("sjmisc")
#install.packages("broom")
#install.packages("haven")
#install.packages("png")
#install.packages("grid")
#install.packages("srvyr")
#install.packages("ggpubr")

library(nnet)
library(effects)
library(magrittr)
library(ggplot2)
library(tidyverse)
library(sjmisc)
library(broom)
library(haven)
library(png)
library(grid)
library(srvyr)
library(ggpubr)

## OPEN DATASET

ca2019.data <- sjlabelled::read_stata("ca2019_ces_poll_regulation.dta")

## DEFINE VARIABLES AS FACTOR OR NUMERIC

ca2019.data$ban3 <- as.factor(ca2019.data$ban3)
ca2019.data$partisan_ahead <- as.factor(ca2019.data$partisan_ahead)
ca2019.data$confusing <- as.factor(ca2019.data$confusing)
ca2019.data$female <- as.factor(ca2019.data$female)
ca2019.data$franco <- as.factor(ca2019.data$franco)
ca2019.data$province <- as.factor(ca2019.data$province)
ca2019.data$postsecondary <- as.factor(ca2019.data$postsecondary)
ca2019.data$interest_high <- as.factor(ca2019.data$interest_high)
ca2019.data$pid_strength <- as.factor(ca2019.data$pid_strength)
ca2019.data$pid_strong <- as.factor(ca2019.data$pid_strong)

ca2019.data$interest <- as.numeric(ca2019.data$interest)
ca2019.data$ideology <- as.numeric(ca2019.data$ideology)
ca2019.data$z_age <- as.numeric(ca2019.data$z_age)
ca2019.data$education <- as.numeric(ca2019.data$education)

## FIT MODEL

ci.fnc = function(effects.obj) {
  col = sym(names(effects.obj$x)) # Capture the name of the effects term
  cis = with(effects.obj, cbind(x, prob, upper.prob, lower.prob)) %>% #Generate a tidy data frame with predictions, lower and upper confidence intervals
    pivot_longer(cols=-col) %>% 
    separate(name, into=c('type', 'response.level'), sep="prob.X") %>% 
    mutate(response.level=gsub("\\.", " ", response.level)) %>% 
    mutate(type=gsub("\\.", "", type)) %>% 
    spread(type, value)
  return(cis)
}

fit.main <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + 
z_age + education + interest + ideology + province, ca2019.data)
summary(fit.main)
nrow(residuals(fit.main)) #Get number of observations

pred.main <- effect("partisan_ahead", fit.main) #Get predicted probabilities

supp.labs <- c("Disagree", "Don’t know", "Agree")
names(supp.labs) <- c("1", "2", "3")

tiff("pred_main.tiff", units="in", width=6, height=4, res=300)

pred.main.graph = ci.fnc(pred.main) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=9),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,.6), breaks = seq(0,.6,.1), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.main.graph

dev.off()

## FIT MODEL BY PARTY

# Liberal Party 

fit.lpc <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
education + interest + ideology + province, data=subset(ca2019.data,pid==1))
summary(fit.lpc)
nrow(residuals(fit.lpc)) #Get number of observations

pred.lpc <- effect("partisan_ahead", fit.lpc) #Get predicted probabilities

tiff("pred_lpc.tiff", units="in", width=6, height=4, res=300)

pred.lpc.graph = ci.fnc(pred.lpc) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.lpc.graph

dev.off()

# Conservative Party

fit.cpc <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
education + interest + ideology + province, data=subset(ca2019.data,pid==2))
summary(fit.cpc)
nrow(residuals(fit.cpc)) #Get number of observations

pred.cpc <- effect("partisan_ahead", fit.cpc) #Get predicted probabilities

tiff("pred_cpc.tiff", units="in", width=6, height=4, res=300)

pred.cpc.graph = ci.fnc(pred.cpc) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.cpc.graph

dev.off()

# New Democratic Party

fit.ndp <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
education + interest + ideology + province, data=subset(ca2019.data,pid==3))
summary(fit.ndp)
nrow(residuals(fit.ndp)) #Get number of observations

pred.ndp <- effect("partisan_ahead", fit.ndp) #Get predicted probabilities

tiff("pred_ndp.tiff", units="in", width=6, height=4, res=300)

pred.ndp.graph = ci.fnc(pred.ndp) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.ndp.graph

dev.off()

## FIT MODEL BY EDUCATION

# No postsecondary education

fit.nopost <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
interest + ideology + province, data=subset(ca2019.data,postsecondary==0))
summary(fit.nopost)
nrow(residuals(fit.nopost)) #Get number of observations

pred.nopost <- effect("partisan_ahead", fit.nopost) #Get predicted probabilities

tiff("pred_nopost.tiff", units="in", width=6, height=4, res=300)

pred.nopost.graph = ci.fnc(pred.nopost) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.nopost.graph

dev.off()

# Postsecondary education

fit.post <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
interest + ideology + province, data=subset(ca2019.data,postsecondary==1))
summary(fit.post)
nrow(residuals(fit.post)) #Get number of observations

pred.post <- effect("partisan_ahead", fit.post) #Get predicted probabilities

tiff("pred_post.tiff", units="in", width=6, height=4, res=300)

pred.post.graph = ci.fnc(pred.post) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.post.graph

dev.off()

## FIT MODEL BY POLITICAL INTEREST

# Low political interest

fit.low <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
ideology + education + province, data=subset(ca2019.data,interest_high==0))
summary(fit.low)
nrow(residuals(fit.low)) #Get number of observations

pred.low <- effect("partisan_ahead", fit.low) #Get predicted probabilities

tiff("pred_low.tiff", units="in", width=6, height=4, res=300)

pred.low.graph = ci.fnc(pred.low) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.low.graph

dev.off()

# High political interest

fit.high <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
ideology + education + province, data=subset(ca2019.data,interest_high==1))
summary(fit.high)
nrow(residuals(fit.high)) #Get number of observations

pred.high <- effect("partisan_ahead", fit.high) #Get predicted probabilities

tiff("pred_high.tiff", units="in", width=6, height=4, res=300)

pred.high.graph = ci.fnc(pred.high) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.high.graph

dev.off()

## FIT MODEL BY PID STRENGTH

# Not very/fairly strongly

fit.notstrong <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
education + interest + ideology + province, data=subset(ca2019.data,pid_strong==0))
summary(fit.notstrong)
nrow(residuals(fit.notstrong)) #Get number of observations

pred.notstrong <- effect("partisan_ahead", fit.notstrong) #Get predicted probabilities

tiff("pred_notstrong.tiff", units="in", width=6, height=4, res=300)

pred.notstrong.graph = ci.fnc(pred.notstrong) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.notstrong.graph

dev.off()

# Very strongly

fit.strong <- multinom(ban3 ~ partisan_ahead + confusing + female + franco + z_age + 
education + interest + ideology + province, data=subset(ca2019.data,pid_strong==1))
summary(fit.strong)
nrow(residuals(fit.strong)) #Get number of observations

pred.strong <- effect("partisan_ahead", fit.strong) #Get predicted probabilities

tiff("pred_strong.tiff", units="in", width=6, height=4, res=300)

pred.strong.graph = ci.fnc(pred.strong) %>% 
  ggplot(aes(partisan_ahead, V1)) +
  geom_pointrange(aes(ymin=L, ymax=U), fatten=1, fill="black", alpha=1) +
  facet_grid(. ~ response.level, labeller = labeller(response.level = supp.labs)) +
  theme_bw() +
  theme(plot.title=element_text(size=10,hjust=0.5),
        axis.title=element_text(size=9),
        axis.text=element_text(size=8),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor = element_line(size = .3), panel.grid.major = element_line(size = .3),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.2), expand=c(0,0)) +
  scale_x_discrete(labels=c("1" = "No", "2" = "None/DK", "3" = "Yes")) +
  labs(y="Predicted Probability", x="Does PID Matches Perceived Party Ahead?", title="")
pred.strong.graph

dev.off()

## COEFFICIENT PLOT

tidy.fit.coef <- tidy(fit.main, conf.int = TRUE, conf.level = 0.95)
tidy.fit.coef

tidy.fit.rrr <- tidy(fit.main, exponentiate = TRUE)
tidy.fit.rrr

tidy.fit <- merge(tidy.fit.coef,tidy.fit.rrr,by=c("y.level","term"))

tidy.fit$estimate.y <- sprintf("%.2f", round(tidy.fit$estimate.y, 2)) 

write.csv(tidy.fit, "tidy_fit_main.csv")

tiff("coef.tiff", units="in", width=6, height=4, res=300)

coef.data <- read_csv("tidy_fit_main.csv")

coef.data <- coef.data [(!(coef.data$term=="(Intercept)") & !(coef.data$term=="province2")) 
& !(coef.data$term=="province3") & !(coef.data$term=="province4")
& !(coef.data$term=="province5") & !(coef.data$term=="province6")
& !(coef.data$term=="province7") & !(coef.data$term=="province8")
& !(coef.data$term=="province9") & !(coef.data$term=="province10"),]

coef.plot <- coef.data %>%
  mutate(y.level = case_when(
    y.level == "2" ~ "Don’t know",
    y.level == "3" ~ "Agree"),
    term = case_when(
      term == "partisan_ahead1" ~ "PID=Ahead? No",
      term == "partisan_ahead2" ~  "PID=Ahead? DK",
      term == "female1" ~ "Female",
      term == "franco1" ~ "Francophone",
      term == "z_age" ~ "Age",
      term == "education" ~ "Education",
      term == "interest" ~ "Political interest",
      term == "ideology" ~ "Left–Right Ideology",
      term == "confusing" ~ "Politics confusing"))%>%
  ggplot(aes(y = term, x = estimate.x, pch = y.level,label=sprintf("%.2f", round(estimate.y,2)))) +
  geom_point(aes(y = term, x=estimate.x), color= "#FF6666", show.legend = FALSE) +  
  geom_errorbarh(aes(xmax = conf.high, xmin = conf.low, height = .12), color ="#FF6666",size = 0.5) +  
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_shape_manual(values = c(16,16)) +
  geom_text(size = 2, nudge_x = .35,vjust = -1) + 
  facet_grid(.~y.level) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) +
  scale_y_discrete(name ="Covariates") +
  scale_x_continuous(name ="Regression Coefficients with Relative Risk Ratios", limits=c(-2,2)) +
  theme(legend.position = "bottom")
coef.plot

dev.off()

## FIT MODEL WITH INTERACTION

ca2019.data = ca2019.data %>% mutate(partisan_ahead = relevel(partisan_ahead, 3))
ca2019.data <- to_factor(ca2019.data, partisan_ahead, postsecondary, interest_high, pid_strong)

ca2019.data$partisan_ahead <- ordered(ca2019.data$partisan_ahead,
levels = c(1,2,3),
labels = c("No", "None/DK", "Yes"))

ca2019.data$ban3 <- ordered(ca2019.data$ban3,
levels = c(1,2,3),
labels = c("Disagree", "DK", "Agree"))

ca2019.data$postsecondary <- ordered(ca2019.data$postsecondary,
levels = c(0,1),
labels = c("No", "Yes"))

ca2019.data$interest_high <- ordered(ca2019.data$interest_high,
levels = c(0,1),
labels = c("Low", "High"))

ca2019.data$pid_strong <- ordered(ca2019.data$pid_strong,
levels = c(0,1),
labels = c("Not Very/Fairly", "Very"))

ca2019.data <- ca2019.data %>% rename(Ban = ban3)
ca2019.data <- ca2019.data %>% rename(Strength = pid_strong)
ca2019.data <- ca2019.data %>% rename(Postsecondary = postsecondary)
ca2019.data <- ca2019.data %>% rename(Interest = interest_high)

# PID strength

int.strength <- multinom(Ban ~ partisan_ahead*Strength + confusing + female + franco + z_age + 
education + interest + ideology + province, ca2019.data, Hess = TRUE)
summary(int.strength)

tiff("int_strength.tiff", units="in", width=6, height=5, res=300)

plot(allEffects(int.strength), colors = "black", selection=9, main="Partisan Effect × PID Strength",
ylab="Predicted Probability", xlab="Does PID Matches Perceived Party Ahead?")

dev.off()

# Postsecondary education

int.post <- multinom(Ban ~ partisan_ahead*Postsecondary + confusing + female + franco + z_age + 
interest + ideology + province, ca2019.data, Hess = TRUE)
summary(int.post)

tiff("int_post.tiff", units="in", width=6, height=5, res=300)

plot(allEffects(int.post), colors = "black", selection=8, main="Partisan Effect × Postsecondary Education",
ylab="Predicted Probability", xlab="Does PID Matches Perceived Party Ahead?")

dev.off()

# Political interest

int.interest <- multinom(Ban ~ partisan_ahead*Interest + confusing + female + franco + z_age + 
education + ideology + province, ca2019.data, Hess = TRUE)
summary(int.interest)

tiff("int_interest.tiff", units="in", width=6, height=5, res=300)

plot(allEffects(int.interest), colors = "black", selection=8, main="Partisan Effect × Political Interest",
ylab="Predicted Probability", xlab="Does PID Matches Perceived Party Ahead?")

dev.off()

## 2019 VOTE INTENTIONS

ca2019.polls <- read_dta("ca2019_campaign_polls.dta")

# Define party colours

colors <- c("lib" = "#d92026", "con" = "#0f2d52", "new" = "#f58320", 
"bloc" = "#00a5ec", "grn" = "#3d953b", "pop" = "#442d7b")

# Graph vote intentions

tiff("ca2019_polls.tiff", units="in", width=6, height=5, res=300)

ca2019.polls.graph <- ggplot(ca2019.polls, aes(x = distance)) + 
geom_point(aes(y = lpc), color = "#d92026", alpha = 1/4) + 
geom_point(aes(y = cpc), color = "#0f2d52", alpha = 1/4) +
geom_point(aes(y = ndp), color = "#f58320", alpha = 1/4) +
geom_point(aes(y = bq), color = "#00a5ec", alpha = 1/4) +
geom_point(aes(y = gpc), color = "#3d953b", alpha = 1/4) +
geom_point(aes(y = ppc), color = "#442d7b", alpha = 1/4)
ca2019.polls.graph

ca2019.polls.graph <- ca2019.polls.graph + geom_smooth(aes(y = lpc, color = "lib"), se = FALSE, method = loess, span = 0.1, size = 1) + 
geom_smooth(aes(y = cpc, color = "con"), se = FALSE, method = loess, span = 0.1, size = 1) + 
geom_smooth(aes(y = ndp, color = "new"), se = FALSE, method = loess, span = 0.1, size = 1) + 
geom_smooth(aes(y = bq, color = "bloc"), se = FALSE, method = loess, span = 0.1, size = 1) + 
geom_smooth(aes(y = gpc, color = "grn"), se = FALSE, method = loess, span = 0.1, size = 1) +
geom_smooth(aes(y = ppc, color = "pop"), se = FALSE, method = loess, span = 0.1, size = 1) +
labs(color = "") + scale_color_manual(labels = c("BQ", "CPC", "GPC", "LPC", "NDP", "PPC"), 
values = colors) + guides(color=guide_legend(nrow=2, byrow=TRUE)) + theme(legend.position = c(.97, 1.02), 
legend.justification = c("right", "top"), legend.box.just = "right", legend.margin = margin(0, 0, 0, 0),
legend.background = element_rect(fill="NA"),legend.text=element_text(size=7), legend.key.size = unit(.4, 'cm'))
ca2019.polls.graph

ca2019.polls.graph <- ca2019.polls.graph + ylab("") + xlab("Number of Days Before the Election") + 
ggtitle("(a) Voting Intention, 2019 Canadian Federal Election") + 
theme(panel.background = element_rect(fill = "white", colour = "white"),
plot.title = element_text(size = 12, hjust = 0.5), panel.grid.major = element_blank(), 
panel.grid.minor = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1),
axis.title = element_text(size = 12), axis.text = element_text(size = 11),
axis.ticks.x=element_line(linewidth=.5),
axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0))) + 
scale_x_reverse("Number of Days Before the Election", limits=c(40,0), breaks = seq(0,40,5), expand = expansion(mult = c(0.02, 0.03))) +
scale_y_continuous("Percentage", limits=c(0,50), breaks = seq(0,50,10), expand = expansion(mult = c(0.03, 0.03))) +
geom_point(aes(x = 0, y = 33.12), colour = "#d92026", size = 2) +
geom_point(aes(x = 0, y = 33.12), shape = 5, colour = "#d92026", size = 4, stroke = 1.) +
geom_point(aes(x = 0, y = 34.34), colour = "#0f2d52", size = 2) +
geom_point(aes(x = 0, y = 34.34), shape = 5, colour = "#0f2d52", size = 4, stroke = 1.) +
geom_point(aes(x = 0, y = 15.98), colour = "#f58320", size = 2) +
geom_point(aes(x = 0, y = 15.98), shape = 5, colour = "#f58320", size = 4, stroke = 1.) +
geom_point(aes(x = 0, y = 7.63), colour = "#00a5ec", size = 2) +
geom_point(aes(x = 0, y = 7.63), shape = 5, colour = "#00a5ec", size = 4, stroke = 1.) +
geom_point(aes(x = 0, y = 6.55), colour = "#3d953b", size = 2) +
geom_point(aes(x = 0, y = 6.55), shape = 5, colour = "#3d953b", size = 4, stroke = 1.) +
geom_point(aes(x = 0, y = 1.62), colour = "#442d7b", size = 2) +
geom_point(aes(x = 0, y = 1.62), shape = 5, colour = "#442d7b", size = 4, stroke = 1.) +
geom_segment(aes(x = 23, y = 0, xend = 23, yend = 48), linetype=2) +
geom_segment(aes(x = 31, y = 0, xend = 31, yend = 48), linetype=2) +
geom_text(size=3,aes(x = 23, y = 49.2, label = "28 Sep")) +
geom_text(size=3,aes(x = 31, y = 49.2, label = "20 Sep"))
ca2019.polls.graph

dev.off()

## MISMATCH BETWEEN PID AND BELIEF ABOUT LEADING PARTY

# Load party logos

lpc <- readPNG("C:/Users/Utilisateur/Desktop/ijpor/lpc.png")
lpc <- rasterGrob(lpc, interpolate=TRUE)

cpc <- readPNG("C:/Users/Utilisateur/Desktop/ijpor/cpc.png")
cpc <- rasterGrob(cpc, interpolate=TRUE)

ndp <- readPNG("C:/Users/Utilisateur/Desktop/ijpor/ndp.png")
ndp <- rasterGrob(ndp, interpolate=TRUE)

bq <- readPNG("C:/Users/Utilisateur/Desktop/ijpor/bq.png")
bq <- rasterGrob(bq)

gpc <- readPNG("C:/Users/Utilisateur/Desktop/ijpor/gpc.png")
gpc <- rasterGrob(gpc)

# Graph mismatch

mismatch.data <- sjlabelled::read_stata("ca2019_ces_poll_regulation.dta")

mismatch.data <- mismatch.data %>% drop_na(pid,partisan_ahead) %>% as_survey() %>% 
  group_by(partisan_ahead,pid) %>% summarize(n = survey_total())

mismatch.data <- mismatch.data %>% group_by(pid) %>% mutate(percent.mismatch = n/sum(n))

mismatch.data <- mismatch.data %>% filter(partisan_ahead==1)
mismatch.data <- mismatch.data %>% filter(pid==1 | pid==2 | pid==3 | pid==4 | pid==5)

mismatch.data <- ungroup(mismatch.data)

mismatch.data$pid <- as.factor(mismatch.data$pid)
mismatch.data$partisan_ahead <- as.factor(mismatch.data$partisan_ahead)

mismatch.data <- mismatch.data %>% add_row(partisan_ahead = "0", pid = "6", n = 1485, 
n_se = 0, percent.mismatch = 0.41783905458)

tiff("mismatch.tiff", units="in", width=6, height=5, res=300)

mismatch.graph <- ggplot(mismatch.data, aes(x = pid, y = percent.mismatch*100)) + labs(x = "Party Identification", y = "") +
geom_col(show.legend = FALSE, aes(color = pid, fill = pid), position = position_dodge(0.4), width = 0.7) +
scale_color_manual(name="Electoral Status", labels=c("LPC","CPC","NDP","BQ","GPC", "All"),values = c('#d92026', '#0f2d52', '#f58320', '#00a5ec', '#3d953b', '#999999')) + theme_bw() +
scale_fill_manual(name="Electoral Status", labels=c("LPC","CPC","NDP","BQ","GPC", "All"),values = c('#d92026', '#0f2d52', '#f58320', '#00a5ec', '#3d953b', '#999999')) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),legend.key.size = unit(.5, "cm"),axis.title.y = element_text(margin = margin(t = 0, r = -4, b = 0, l = 0)),legend.text = element_text(size=10), 
legend.title = element_text(size = 10),plot.title = element_text(size = 12, hjust = 0.5), 
axis.title = element_text(size = 12), axis.text = element_text(size = 11)) + 
scale_x_discrete(labels=c("LPC","CPC","NDP","BQ","GPC", "All")) + scale_y_continuous(limits=c(0, 80), 
breaks = seq(0, 80, 20)) +
ggtitle("(b) Mismatch Between PID and Perception of Leading Party") + 
theme(axis.title.x = element_text(margin = margin(t = 8, r = 0, b = 0, l = 0)), axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),plot.title = element_text(size = 12, hjust = 0.5)) +
geom_text(size=3,aes(x = 1, y = 24, label = "(N = 351)")) +
geom_text(size=3,aes(x = 2, y = 24.5, label = "(N = 312)")) +
geom_text(size=3,aes(x = 3, y = 73, label = "(N = 482)")) +
geom_text(size=3,aes(x = 4, y = 71, label = "(N = 156)")) +
geom_text(size=3,aes(x = 5, y = 69, label = "(N = 184)")) +
geom_text(size=3,aes(x = 6, y = 45, label = "(N = 1485)"))
mismatch.graph

mismatch.graph <- mismatch.graph + annotation_custom(lpc, xmin=0.93, xmax=1.12, ymin=1, ymax=31)
mismatch.graph <- mismatch.graph + annotation_custom(cpc, xmin=1.9, xmax=2.1, ymin=1, ymax=31)
mismatch.graph <- mismatch.graph + annotation_custom(ndp, xmin=2.75, xmax=3.25, ymin=1, ymax=129)
mismatch.graph <- mismatch.graph + annotation_custom(bq, xmin=3.87, xmax=4.12, ymin=1, ymax=122)
mismatch.graph <- mismatch.graph + annotation_custom(gpc, xmin=4.88, xmax=5.14, ymin=1, ymax=119)
mismatch.graph

dev.off()

## COMBINE VOTE INTENTIONS AND MISMATCH GRAPHS

polls.mismatch <- ggarrange(ca2019.polls.graph, mismatch.graph, nrow=1, ncol=2) 

ggsave(plot = polls.mismatch,filename="polls_mismatch.png", units="in", width=10, height=5, dpi=300)